##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(ArchR)
library(BSgenome.Hsapiens.UCSC.hg38)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--cluster_type_file"), type = "character"),
    make_option(c("--cluster"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    comine_data_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined_peak.combineRNA.qc.Rdata"

    ## 细胞类型注释文件
    cluster_type_file <- "~/20231121_singleMuti/config/cluster_celltype.csv"

    cluster <- "cluster0"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/subcell"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
cluster_type_file <- opt$cluster_type_file
cluster <- opt$cluster
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

a <- load(comine_data_file)
atac_proj <- testis_combined_peak_combineRNA

## 分细胞类型跑
cluster_type <- fread(cluster_type_file , header = F)
celltype_name <- subset( cluster_type , V1 == cluster )$V2

###########################################################################################
## 提取细胞类型
subCells <- getCellNames(testis_combined_peak_combineRNA)[as.character(testis_combined_peak_combineRNA@cellColData[["cell_type"]]) %in% celltype_name]

## 提取特定细胞的
sub_proj <- subsetArchRProject(
  ArchRProj = testis_combined_peak_combineRNA,
  cells = subCells,
  outputDirectory = out_path,
  dropCells = TRUE,
  force = TRUE
)

###########################################################################################
## 基于chromvar鉴定motif
projHeme5 <- sub_proj

projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "cisbp", name = "Motif" , force = TRUE)

projHeme5 <- addBgdPeaks(projHeme5 , force = TRUE)

projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "Motif",
  force = TRUE
)

###########################################################################################
## scATAC-seq 和 scRNA-seq 数据识别峰到基因链接
## projHeme5 <- addPeak2GeneLinks(ArchRProj = projHeme5,useMatrix = "GeneExpressionMatrix")
atac_proj <- projHeme5

###########################################################################################
out_file <- paste0( out_path , "/" , cluster , ".combineRNA.motif_peak2gene.Rdata" ) 
save( atac_proj, file = out_file )
saveArchRProject( atac_proj , out_path )
